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Neutron star crustal matter, whose properties are relevant in many models aimed at explaining 
observed astrophysical phenomena, has so far always been studied using a mean field approach. In 
order to check the results obtained in this way, a sensible next step is to make use of a realistic 
nuclear potential. The present paper extends the periodic-box Fermi HyperNetted Chain method to 
include longitudinal-isospin dependence of the correlations, making feasible a study of asymmetric 
crustal matter. Results are presented for the symmetry energy, the low-density neutron star equation 
of state and the single particle neutron and proton energies. 

PACS numbers: 



I. INTRODUCTION 

Various observed phenomena connected with neutron stars (NS) - such as glitches - are thought to be largely 
determined by their crustal properties. Glitches are abrupt changes in pulsar periods, [3311 widely thought to be 
related either to starquakes or to vortex pinning, although there are also other possibilities A key role in these 
interruptions in the usually regular behaviour of the pulsar spin period is thought to be played by the interaction 
between superfiuid neutrons and normal matter in the crust. 

Although a number of recent papers have focused on studying NS crustal properties (e.g. [| d), none of them have 
gone beyond a mean field approach for describing the nuclear interaction. That sort of approach allows one to give 
a self-consistent description of all of the nucleons forming the crustal matter, either as part of the large neutronized 
neutron-drip nuclei which form the lattice or as part of the neutron superfiuid which flows through it. (For a recent 
review of NS crustal matter see However, major advances have been made in recent years in developing nuclear 
many-body methods for dealing with strongly interacting fermions, which have enabled very accurate calculations to 
be made of several properties of nucleonic matter at low temperature, in both normal and superfiuid phases, which 
have clearly shown that NN correlations play a fundamental role and cannot be disregarded. The short-range behavior 
of these arises from the strong repulsion of the nuclear potential at short inter-nucleon distances, and the long-range 
behavior results from the tensor interaction coming from pion-exchange. These types of behavior have proved to be 
a distinctive feature of NN correlations and hence of the nuclear medium, and they lead to potentially measurable 
effects related to NS structure and the neutrino mean free pathQ. 

These technical advances make it possible to perform ab initio calculations of the structural properties of the NS 
nuclear medium, all the way from the crust to the inner core, fully based on a realistic bare NN interaction. This 
paper is concerned with addressing this challenging problem. 

Obviously, ab initio calculations have a much more limited range of applications than calculations with mean field 
theories which, however, are based on ad hoc effective interactions, which are supposed to include the main features 
of the NN correlations in some average way. Therefore, a second important goal of these studies is to derive mean 
field effective interactions from first principles, starting from a unique bare nuclear interaction. 

We base our investigation here on Orthogonal Correlated Basis theory (OCB) and nuclear Quantum Monte 

Carlo (QMC) methods, particularly Variational Monte Carlo (VMC) and Auxiliary Field Diffusion Monte Carlo 
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(AFDMC)0. All of these methods enable one to make very accurate calculations of nuclear matter interacting with 
modern NN potentials, but they have all only been used so far for homogeneous matter, for nuclei or for neutron 
droplets. The NS crustal matter is composed of a lattice of neutron-rich nuclei surrounded by superfluid neutrons. 
Typical mixtures are characterized by a lattice length of ~ 20 — 60 fm and density values of the neutron soup ranging 
from ~ 0.001 fm~ 3 to ~ 0.1 fm~ 3 . In spite of having such relatively low densities, the neutrons are still in a regime of 
strongly-correlated fermions, because the neutron-neutron scattering length is large, a nn ~ — 20 fm and the resulting 
values of kpa nn range from ~ —6 to ~ —30. 

As far as QMC is concerned, one has to deal with periodic boxes of dimension ~ 20 — 60/m, with an heavy nucleus 
at the center and a few thousand neutrons surrounding it. This requires very resource-intensive simulations, which 
cannot be massively-parallel. Because of this, OCB theory will be used first in order to get optimal variational wave 
functions, and variational estimates of the binding energy per particle (as a function of the the neutron density and of 
the symmetry parameter a = (N — Z)/A, as well as for computing quantities like energy spectra, response functions, 
spectral functions, etc., which are at present beyond the reach of QMC simulations. 

In pursuing these goals, we are faced with two kinds of problems: (i) the lattice structure of the NS matter does 
not allow for calculations in the thermodynamic limit: systems with 2 — 3 thousand neutrons are still far away from 
that; (ii) we are dealing with asymmetric matter having N > Z and we cannot just rely on using the two extreme 
cases a = and a — 1 and then making a quadratic interpolation between them for all of the other cases, because of 
the presence of the nucleus in the box. 

The first of these problems requires relying on the existing PB-FHNC version of the FHNC theory. The second one 
requires us to rewrite the PB-FHNC to allow an isospin dependence. 

The present paper is devoted to clarifying these two points, which represent fundamental methodological steps 
towards making a truly microscopic and unified treatment of NS crustal matter. For doing this, we have extended the 
PB-FHNC method to deal with homogeneous asymmetric nuclear matter described by correlated basis functions, in 
which the correlation operators have a longitudinal isospin dependence. Further considerations limit the variational 
choice of the isospin dependence of the correlation operator F to its longitudinal component t z (i)t z (J) only. It is then 
important to ascertain how good is such a variational choice compared with the full one given by r(i) • t(J). As a 
first application, we have considered simple two-body nuclear potentials of the v± type, which do not include tensor 
components but have full spin-isospin dependence. Calculations have been made of the equation of state (EOS) and 
of the single particle spectra at various values of symmetry parameter. 

The results obtained are very encouraging. The iterative process developed for solving the new PB-FHNC equations 
converges rapidly and gives stable solutions. Moreover, the longitudinal isospin dependence is able to account for more 
than 80% of the full isospin dependence for the interaction models which we have considered, irrespective of the values 
of the density p and the symmetry parameter a. Interesting results are obtained for the symmetry energy, particularly 
in the low-density region. We compare these results with the Bethe-Brueckner-Goldstone (BBG) calculations of ref. 

M- 

Since we want to address this paper and subsequent related ones to the astrophysical community, we repeat here 
some material which has already been published in the nuclear physics literature, in order to make the presentation 
comprehensible. We recognise that readers who are not specialists in nuclear physics will need to be strongly motivated 
in order to work through all of this, but we are aiming here to provide a "bridge" for those strongly-motivated readers. 

The paper is organized as follows. In the next sections we will discuss the nuclear interaction, focusing on the two 
particular interactions that have been used in our calculations. Then, in section ITTT1 we describe the state-dependent 
particle box FHNC scheme used in this work, and finally we present our results and give conclusions. 



II. NUCLEAR POTENTIAL 

A realistic nuclear potential is usually written as a two-body potential (e.g. Argonne V18 (ll|) plus a three-body 
contribution (e.g. Urbana IX [IH) which becomes increasingly important beyond half of the nuclear saturation density 
(po = 0.16 fm" 3 = 2.7 x 10 14 gcm" 3 ) (see [jj] and references therein). It has been shown that in medium n-body 
potentials, with n > 2, can be successfully simulated by two-body density dependent terms (i~3l. Hi}. 

In this paper, we consider the Illinois class of two-body potentials, which are characterized by having a strong local 
contribution, given by their first six spin-isospin-dependent components: 

v 6 (iJ) = J2v (p) (r ij )d%\iJ), (1) 
P =i 

with: 

Off 1 - 6 = (l,< Hj ,S ij )®(l,T ij ), (2) 
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FIG. 1: AT4' potential components 



where Ty = r,; • Tj and <7jj = cti ■ <jj , with r and a being the the Pauli matrices acting respectively on the isospin and 
the spin of a nucleon, and Sy = (3f^f^- — 5 a p)<j™<j? is the tensor operator. 

The most important non-local components of the Illinois potentials are the spin-orbit terms, Ljj -Sy and -S^Ty, 
where Ly and are the relative angular momentum and the total spin of the nucleon pair. Model potentials that 
include the spin-orbit components are denoted as v%. Other components include L 2 , (Ly • Sy) and symmetry 
breaking operators, giving the total of 18 components of Argonne Vi$ (AV18) . 

In its full form, AV18 gives an almost perfect fit to the NN data up to the meson production threshold. Other 
realistic NN potentials, such as the Bonn and Nijmegen potentials, fit the NN data equally well. All of these realistic 
interactions are more or less equivalent in describing the properties of light nuclei but the latter ones are basically non- 
local, and so they are much more difficult to handle with many-body theories such as OCB and AFDMC. Moreover, 
NN correlations are much better described in r-space, in which they are clearly distinguished from relativistic effects. 

There is strong evidence that the first 8 components of the Illinois-type potentials are sufficient for giving a realistic 
description of the nuclear medium. Such model potentials can be obtained by simply truncating AV18 after the first 
8 components. A better choice, however, is to keep the vs form and re-fit the NN data. That has been done (llj . 
and the corresponding model interaction is known in the literature as Argonne Ug/ or AV8'. The differences between 
AV18 and AV8' are quite small and can safely be treated perturbatively. 

Other widely used interaction models from the same class of potentials are v& , which has the form of cq. [T] and v^> , 
in which the tensor components are also omitted. The AV6' and AV4' potentials can be found in ref.flG]. They should 
really be considered only as toy potentials, although they can reproduce a certain amount of nuclear physics data. 
They are certainly very useful, though, for checking many-body techniques and for finding the relative importance of 
the tensor and spin-orbit correlations. 

One of the potentials which we have used in our calculations is the S3 potential derived by Afnan and Tang by fitting 
low energy nucleon-nucleon s-wave scattering data (up to 60 MeV) [T3|. This is a potential of the Serber type, 
and is therefore defined for only the even states. It provides a reasonable description of some basic properties of H 3 
and He 4 such as, for instance, the binding energy and the root-mean-square radii. As in the PB-FHNC calculations 
of ref. [l8[, we have added to the original S3, an interaction for the odd channels given by the repulsive part of the 
even channels. The four components of the resulting potential, which we denote as the AT4' potential, are given by 

v {1) = v c = v R + — (v A t(rij) + v As (rij)) , 

f (2) = v T = -—(3v A t{rij) -VAs{nj)) , 

u (3) = v<j = jx{vAt(rij) - 3v A s(n 3 )) , 

u (4) = v aT = -—(v A t{rij) + VA S (rij)) , (3) 
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with the repulsive and attractive components given by 

Vr = w s iexp(-^ir 2 ), 



3 

VAs = y^u si ;exp(-/3 st r 2 ). 

i=2 
3 



VAt 



= ^u ti exp(-ft 4 r 2 ), (4) 



where the strengths v si and v t i and the gaussian coefficients j3 S i and Pu are given in ref. [l7|. The four components 
of the AT4' potentials are compared in Fig. 1. 



The PB-FHNC calculations of ref. [18[ were performed for the case of Jastrow correlated theory. The resulting 
variational energy for a system of 2060 nucleons at the experimental equilibrium density po is — 15.150MeV, which is 
not too far away from the FHNC/SOC result, E = -16. 184 MeV . 

In order to compare with the Jastrow results of ref. fl8T | , we have also performed variational calculations with the 
potential t>4, obtained by truncating AV8' after its first four components, and denoted here as {AV%')±. 

III. CORRELATED BASIS FUNCTIONS 



The correlated basis functions for a strongly-correlated Fermi fluid are given by 

F\n] 



F^F I n}? ' 



(5) 



where | n] is a generic eigenfunction of the Fermi Gas hamiltonian and F is a correlation operator. The label n 
indicates the number of particle-hole excitations and | 0] denotes the Fermi Gas ground state. The set of Fermi 
Gas states | n] is orthonormal, whereas that of the correlated states | n) is not, because the correlation operator F 
breaks the orthogonality. We need to restore orthogonality by following the two-step procedure outlined in ref. H: 
orthonormal correlated states are denoted as | n) . 

A. Properties of the orthonormalization process 

It can be shown that the orthonormalization procedure has a number of important properties. Formally in OCB 
theory, the Hamiltonian H is written as the sum of an unperturbed term Hq and an interaction term Hi: 

(n | H I m) = (n \ H \ m)S nm = H nn S nm , 

(n | Hi I m) - (n\H\ m)(l - 5 nm ) = H nm . (6) 

The interaction Hamiltonian Hi is simply the non-diagonal part of H. The main property of the orthogonalization 
process of OCB theory is that the diagonal matrix elements are not modified by it, i.e.: 

H nn = (n \ H \ n) = (n \ H | n) + terms of order I/O, . (7) 

This guarantees that the variational estimates (/„ \ H \ I n ) are maintained after orthogonalization. A second 
important property is expressed by the following equation: 

(n\H\n) = (0|i/|0) + ^e„(p l )-^e„(h l ) 

Pi h, 

+ terms of order 1/S7 . (8) 

In this paper we restrict attention to the diagonal matrix elements of the hamiltonian on the ground state | 0) and 
on the one particle-one hole states | ph), given by 

/ dRA [0i . . . <p A ] F] l F jl A [4> a ... 0i] 
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where (j> n = exp(ik n • Y{]Xn{si)I n {Ti) and the orbitals 4>\. . .(f>A correspond to the Fermi sea states in the case of 
the ground state energy (p = h = k F ) and include the excited state orbital p in place of the Fermi sea state orbital 
h in the case of the particle-hole excitation. A is the antisymmetrization operator. Therefore A [4>i . . . 4>a] is a Slater 
Determinant of plane waves. Note that in a Periodic Box treatment, one has a finite number of nucleons and a fixed 
value for the length L = SI 1 / 3 of the box, given by p p + p n = p = A/Q. 

Integration over dR extends to all of the a coordinates (R = ri, T%, . . . , ya) and includes summation over all of the 
spin and isospin variables. 



B. The correlation operator 

Since we are considering a U4 model interaction, we do not need to include tensor correlations in F. In this case, 
the standard choice is given by 



Fa = S 



n hdj) 

i<i=l 




' 4 

E 

P =i 



/ (p) (r, 



(10) 



" (p) 

where the symmetrization is needed because, in general, the operators do not commute with each other. In the 
calculation of the matrix elements of any given operator, it is not known how all of the orderings of the right hand 
side of eq. (|10p can be taken into account. The best known approximation is the so-called FHNC/SOCp^j|, which has 
been shown to give reliable results in a number of nuclear matter calculations. However, this approximation cannot 
be used for asymmetric nuclear matter because of the presence of the operator r(i) ■ t(J) in F4. Because of this, it is 
important to check the variational relevance of the longitudinal isospin-dependent operator Tz (i,j) = T z (i)r z (j) as 
compared with r(i) ■ r(j). The operators Tz (i,j) commute with each other, and so can easily be used in asymmetric 
nuclear matter calculations. To achieve this goal, we consider the following correlation operator: 



Fj L = Yi [fNN(Tij)PNN(iJ)+fpp{rij)Ppp(i,j) 

i<j=l 

+ fNp(nj)PNp(iJ) + fpN(r%j)PpN(i,j)] , 



(11) 



where the projection operators P a b are given by 



PNN{i,j) 

Ppp{iJ) 
PNp(i,j) 
PpN(i,j) 



1 


+ r z (i) 


1 


+ r z {j) 




2 




2 


1 


-T z (i) 


1 


-r z (j) 




2 




2 


1 




1 


-r z (j) 




2 




2 


1 


~r z {i) 


1 


+ r z {j) 




2 




2 



(12) 



The four scalar functions /nn, fpp, Jnp and fpN have to heal smoothly to 1, thus giving an uncorrelated system, 
for rij greater than a certain healing distance d chosen so as to minimize the ground state energy per particle of the 
system (0\H\0)/A. The results discussed in this paper are obtained under the following assumption: 



fNN(nj) = fpp (nj ) = /||(ry), 

/jvp(ry) = fpN(rij) = f±(rij). (13) 

The correlation functions /|| (r^ ) and f±(rij) are obtained by solving a set of second order differential equations [l9| 
whose detailed application to the system at hand is described in appendix [X] They distinguish isospin parallel from 
isospin antiparallel correlations. Such isospin dependence does not completely resolve the difference between T = 1 
and T = channels as in the case of F4 of eq. (flU)) . Nevertheless, as will be shown in this paper, it still provides a 
very good description of the isospin dependence of nuclear correlations induced by a nuclear two-body potential of 
the V4 type. 
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IV. POWER SERIES EXPANSION AND DIAGRAMMATIC RULES 



The CBF matrix elements given in eq. ([9]) are calculated by first applying standard Fantoni-Rosati (FR) cluster 
expansion techniques [20I l2l| to (ph|i?|ph) and by summing up the resulting series of cluster terms using the FHNC 
integral equation methods. The FR expansion is based on expanding both the numerator and the denominator of 
eq.© in powers of the functions h a (r) given by: 



MO = /aW - 1 . (" = NN, PP, NP, PN) . 



(14) 



The quantity Fj L Fjl appearing in the denominator of eq. © is then decomposed into a series of cluster operators 



as follows: 



F\ T F JL = ^4_Hf (1,2) 



JL r JL 



l+^^ 3 )(l,2;z) 



/ 2 (1,2) = £/»(r ia )P (l,2) 



(15) 



where each cluster operator / 2 (l,2)AfW(l,2;i) is expressed in terms of products of h a functions and t q projection 
operators and correlates the two interacting particles 1 and 2 with the other n — 2 particles in the medium. A 
similar expression is obtained for the quantity P\ l HFjl appearing in the numerator of eq. ([9]), where / 2 (1,2) is 
substituted by /(l, 2)H{\, 2)/(l, 2), with if (1,2) including the potential energy operator 1)4(1,2), and the kinetic 
energy operators: V 2 /(1,2), Vi/(1,2) • Vi/(1,»), Vi/(1,2) • Vi\n], etc. 

Inserting the cluster decomposition of eq. (fl"5)) into eq. ©, each cluster operator / 2 (1, 2)X^ (1, 2; i) gets multiplied 
by the n-body Fermi Gas distribution g n (l, . . . , n). In appendix (|B|), the distribution g n (l, . . . , n) is expressed in terms 
of the uncorrelated one-body density matrix (also called the exchange function) £^(i,j) and (p(i,j) which are given 
by 



^ n— 1 m—up 7 dow7i 

l z 

WW) = -E^W^W^M^) E x™Wx™(j), 



(16) 



m— up, down 



where the sum is extended over the occupied states. They are normalized to pn and pp respectively. It is useful to 
define a four component vector function given by 



Zpp(Tij) = ip(Tij), 

iNp(Tij) = e PN (r lJ )=0. (17) 

The net result is that / 2 (l,2)A?(")(l,2;i) gives rise to a sum of n-body cluster integrals, whose integrands are 
products of dynamical correlation functions ft, Q (ry)and exchange correlations £ a (ry) . 

A very important property of the FR cluster expansion, when applied to finite systems, is that in both the numerator 
and denominator of the diagonal matrix element of eq.©, the summation over the cluster integrals can be extended 
beyond the order A, which is the maximum for a system of A nucleons. In fact it can be extended up to infinity 
because any g n (l, ...,n) built with the exchange functions £ a (rij) given in eq. (|16j) . with n > A vanishes. This 
property enables us to use all of the FR cluster expansion properties which are valid for a system with an unlimited 
number of nucleons, such as for instance nuclear matter. 

More details about the FR decomposition can be found in ref. [22j. We report here only its main properties so as 
to help the reader get a quick understanding of the original papers. 

Each cluster integral is most conveniently represented by a cluster diagram. These diagrams are built by following 
a few convenient rules [2lj |: 

1. each point represents a particle. Filled points - or internal points - represent in-medium particles while unfilled 
points represent external interacting particles. 
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unlinked linked and simple linked and composite 

FIG. 2: Examples of cluster diagrams: A is unlinked, B is nodal, C is elementary and D is composite. All the points correspond 
to neutrons 

2. two points are linked either by a dashed line, representing h a (r) or by a solid oriented line, representing £(r) or 
by both ; 

3. any internal point is reached by one or more dashed lines; two dashed lines cannot be superimposed; 

4. solid lines form closed loops and different loops cannot have any common point. 

5. each internal point carries a proton or neutron density factor, depending on which type of particle it represents; 

6. all of the particles belonging to an exchange loop are in the same spin-isospin state. Each loop (except for those 
comprising two particles) is counted twice, because there is one oriented clockwise and one anti-clockwise. The 
two-particle loops are counted only once. In addition, one has to sum over loops of different spin-isospin states. 
The loop sign is given by (— where n is the number of points on the loop. The global factor C„ of an 
n-particle loop in spin-symmetric matter is given by 

C n = R" +1 x 2 x 2 x (p P + p n N ) , (18) 

where pp and pn are the proton and neutron matter densities respectively. For symmetric nuclear matter 
{fip = Pn = p/2), the factor C n becomes (-) n+1 8(p/4) n . 

There are some exceptions to these rules, which occur in the calculation of the expectation values of the n • r 2 and 
{<j\ ■ o- 2 )(ti • r 2 ) potential terms. These will be discussed later, in connection with the calculation of the energy per 
particle. 

The linked cluster theorem[21] holds also in the case of the JX-correlated basis. This theorem states that non-linked 
cluster diagrams (i.e. diagrams which are built from two or more completely unconnected parts and which diverge in 
the thermodynamic limit) cancel exactly between the numerator and denominator of eq. [9] so that one is then left 
with a series of linked cluster diagrams. 

Examples of cluster diagrams are shown in FigJ5] Diagram (J5^) is unlinked and is therefore forbidden. The remaining 
diagrams are all allowed. In the calculation of the expectation value of the scalar component of the two-body potential 
diagram, (J2}d) corresponds to the following contribution 

Diagram (lb) -> -— p% I dY 1 2dYidY 3 dY k f?AY 1 2)v c (r 12 ) 
^Pn J 

X h^T U )h\\{T 2 ^{T jh ){ ~^^ \ . (19) 

Linked cluster diagrams are subdivided into two classes: simple and composite (or hyper-chain). Simple diagrams 
are further classified as nodal (or chain) or elementary. 

Nodal diagrams (iV-diagrams) are defined as diagrams having one or more nodes: internal (filled) points that are 
necessarily crossed by any path going from one interacting point to the other. For instance, diagram ^jp) is nodal 
and the points labelled with i and j are both nodes. Not nodal diagrams (JT-diagrams) include both the elementary 
(-E-diagrams) and the composite ones. Composite diagrams are obtained by combining two or more nodal diagrams. 
Diagram (J5J1) is composite and is composed of two nodal subdiagrams. Elementary diagrams (-E-diagrams) are 
the remaining ones. They are neither composite nor nodal. They can be constructed by first identifying the basic 
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topological structures we want to include. Such basic structures have the property that each of their internal points 
is reached by at least three links and the two external points by at least two. The links are building blocks given by 
N + X diagrams. The basic structures are characterized by their number of points, the minimum being four (£4). 
Diagram (J2fc) is an example of four-point elementary diagrams. 



V. PB-FHNC SCHEME 

The PB-FHNC integral-equation method[22| gives an easy way of summing the series of cluster terms by specifying 
rules for building diagrams, using other diagrams as building blocks in an iterative way. The method extends Hyper 
Netted Chain (HNC) theory to the case of correlated quantum Fermi systems; HNC theory has been widely used 
in statistical thermodynamics [23j and, more recently, has been applied to CBF calculations for low temperature 
Bose systems such as liquid 4 He and 3 He impurity in 4 He 24]. It is based on two basic algorithms: the chain and 
hyper-chain. 

a. Chain algorithm. This consists of summing up the whole series of TV-diagrams made from a given building 
block. One simply takes the sum of the geometric series of the Fourier Transforms of the building block function, 
which is given by a subset of Jf -functions 

N(k) = pX 2 (k) + p 2 l 3 (k) + • • • = pX( & , (20) 

1 - pX(k) 

which can be expressed in terms of the following integral equation: 



Nfa) = pj dr kj X(r ik )(N(r kj ) + X{r kj )). (21) 

This formula can be interpreted in terms of probabilities. Each integral can be thought of, given two external 
particles as the probability to find a third in-medium particle k - expressed by p - times the probability of i 
interacting with k times the probability of k interacting with j. 

b. Hyper-chain algorithm. This consists of summing up the whole class of X-diagrams made from a given 
subset of iV-diagrams: 



X{ rij ) = /^r^exp^r^+^ry)]-^^)-!, (22) 

where / 2 (ry) is given by exp(—V(jCij)/KT) in statistical thermodynamics and by the correlation function of the 
scalar Jastrow ansatz Yl(f( r ij)) m the variational calculations of zero temperature Bose systems. The function S(ry) 
corresponds to the E-diagrams, which cannot be calculated in a closed form, like the N- and X-diagrams. The 
meaning of X becomes clear if we imagine expanding the exponential in series: we are summing an increasing number 
of N and E diagrams. The result of this procedure will give of course composite and elementary diagrams. In addition 
to those the exponential term includes the nodal diagrams and the identity which we must subtract in order to get 
the sum of the not-nodal diagrams. The function E(rij) is a functional of iV(ry) and X(ry) and in general it is 
approximated by the first few-body basic diagrams (the lowest of which is the four-body basic diagram, like the 
diagrammatic structure underlying diagram^. 

c. FHNC algorithm . Eqs. (|2"Tj) and (|2"2")l arc formally identical to the HNC equations of statistical thermo- 
dynamics and can be solved in an iterative way by means of the following steps: 

1. take N{rij) = 0; 

2. compute X{jCij) using eq. (f2"2"j) and taking the function AT(ry) from the previous step; 

3. compute -/V(r,y ) using eq. (|2"Tj) and taking the functions X(rij) and iV(ry) on the r.h.s. from the steps 2 and 3 
respectively; 

4. return to step 2, and continue until convergence is obtained[34]. 
The pair correlation function 
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FIG. 3: Labeling of external points. See text for further details. 
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FIG. 4: Building of an Ndd- graphical counterpart to eq. 1251 



3(r 12 ) - j^-—^, (23) 



where 



g{r 12 ) = 1 + N{r 12 ) + X{r 12 ) = f 2 (r 12 )exp{N{r 12 ) + E{r 12 j). (24) 

In the case of Fermi systems, in addition to the dynamical correlation bonds, there are also the exchange bonds, 
with the diagrammatic rules given in Section IV. Because of this, the FHNC method requires further subdivision of 
diagrams by labeling the exchange character of the two external points (e.g. fig[3]). Each point in a diagram is labelled 
with a d unless it is reached by an exchange line part of a closed loop (labelled with e). Therefore one has four different 
nodal functions N dd (rij), N de (r i:j ), N ed (r i:j ) and iV ee (ry), with N de (r i:j ) = N ed (rji). Similarly one has four X- and 
E-functions. In addition one needs to introduce another class of functions, those in which the two external points are 
joined by an open loop of exchange lines, which are denoted with the label cc. This last class of points should not 
be present in allowed diagrams; diagrams with c points are non-physical (i.e. they have no physical meaning if taken 
alone) but we include them as useful building blocks, to construct diagram having a closed loop passing through the 
two external points 1 and 2. 

It should be noticed that, as in HNC theory, the external points i and j of N, X and E diagrams summed up at 
a given iteration of the FHNC scheme, may become internal points of the next generation of diagrams. The true 
external points 1 and 2 are those at convergence. 

The chain algorithm for a Fermi system has to take into account the statistical nature of the convolution node, 
which can be either d or e in the case of the chain equations for N dd , N de and N ee , and has to be necessarily of the 
type c for N cc . As an example, the equation used to build N dd (vij) is: 

N dd (r i3 ) = p I dr kj X dd (r lk ) [X dd (r kj ) + N dd (r kj )] 
Jn 

+ p dr kj X dd (r ik ) [X ed (v k] ) + N ed (r kj )} 
Jn 

+ p dr kj X de (r ik ) [X dd (r kj ) + N dd (r kj )] , (25) 
Jn 

where in the first convolution integral the node is of the d type and, in the remaining two, it is of type e . Fig. (|4]) 
exemplifies the above chain equation. 

In close analogy with HNC theory, the equation used to build an X dd diagram is 



X dd (r 12 ) = g dd {r 12 ) - N dd (r 12 ) - 1 , 



(26) 
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where 

9dd{ri 2 ) = f(r 12 ) exp (N dd (r u ) + E dd (r 12 )) . (27) 

The complete set of formulae that are needed to compute the pair distribution function and the energy per particle 
are given in p2l |. At present we do not know any formula able to give us a useful prescription to sum E diagrams 
in a closed form, like those for N and X. Luckily enough their contribution has been shown to be almost negligible 
in the case of translationally invariant nuclear systems [251 ] . This allows us to neglect them for all practical purposes 
and hence to use the so called FHNC/0 approximation. 

A. Vertex corrections 

Linked diagrams can be divided into two classes: (i) reducible diagrams (like diagram b of fig. ^ having one or more 
reducibility points, which are the only contact points of two sub-diagrams; (ii) irreducible diagrams (like diagrams c and 
d of figH]) which have no reducibility points. The cluster integral corresponding to reducible diagrams is factorizable in 
the product of the cluster integrals corresponding to the underlying irreducible diagrammatic structures. The FHNC 
scheme sketched in the previous section sums up the irreducible diagrams only. In fact, they are the only remaining 
diagrams in pure Jastrow theory (/|| = /j_), because all the reducible ones cancel each other exactly. For instance, 
diagrams a and b of fig. [5] cancel each other because the exchange loop insertion leads to a factor —1 due to the 
orthonormality of the single particle orbitals <p n . However the cancellation is no longer true in the more general case of 
the ansatz given in eq. [TU]and in eq. [T3]or when the Slater determinant of plane waves is substituted by a BCS wave 
function to describe a superfluid fcrmi system, like for instance neutron matter at low density and zero temperature. 

It has been proved [21| that one can still use the FHNC scheme previously described, paying the price of renor- 
malizing the various points of the irreducible diagrams by proper vertex corrections, which take into account all the 
possible one-body subdiagrams that can be linked to them. The process is exemplified by the diagrammatic equation 
displayed in fig. [5j There are two types of vertex corrections, £d(r) for points not reached by exchange lines and £ e (r) 
for the others . They are given by the following equations 

&(r) = (l + C/ e (r))exp(t/ d (r)), 

£ e (r) = exp((7 d (r)), (28) 

where U e and U d are the sums of diagrams with e and d starting points respectively. In practice, this is accomplished 
by integrating the functions X dd {rij) and JQ e (ry) over rj for U d and X ed {r^) and X ee (rij) over rj for U e and then 
correcting them to avoid overcounting that would arise because of the increased symmetry passing from a two-body 
diagram to a one-body one. For instance triangular diagrams of fig. [5]) have different symmetry factors. The factor 
is 1 for the two-body diagram ([5ji) and 1/2 for the one-body diagram (JSJs). The derivation of the equations leading 
to U d and U e can be found in ref. [2~fl ] 

In the above equations U d comes in the exponential to account for the unlimited number of U d terms that can stem 
from any point either of type d or type e. On the contrary, one can have one U e term only as a vertex correction of a 
point of type d. It should be noted that £ d corresponds to the sum of all the possible one-body linked diagrams and, 
hence, to the one-body correlation function that, for a homogeneous system, is gi(r) = 1. The sum rule t; d = 1 can 
be used as a measure of the accuracy of the FHNC approximation. 

As a final remark, we should note that vertex corrected (or renormalized) diagrams do not need to obey our third 
FHNC rule: i.e. there may be diagrams with internal points not reached by a dashed line. However this requires 
that the point should have been reached by a dashed line pertaining to a subdiagram accounted for by a correction, 
implying that we will introduce a third correction £ c that includes the same diagrams as £ e except for the identity: 

6 = 6-1, (29) 

After having sketched the main instruments of PB-FHNC theory we extend it in the following to include longitudinal 
isospin dependence in the correlation operator and to the N Z trial functions. 

VI. STATE DEPENDENT PB-FHNC EQUATIONS 

The PB-FHNC equations derived in this paper are obtained for the more general form of the correlation operator, 
given in eq. (fTTj) and not for the restricted one of eq. (fT5)). It follows that nodal, composite, elementary functions, as 
well as distribution functions, will have the structure of a four-component vector 
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a be 




6 b 6 

12 12 

FIG. 5: Diagrammatic exemplification of vertex correction. Diagrams (a-c) are summed up to give £d(ri). Diagrams (d) and 
(e) show a basic property of the FHNC equations for Ud and U e - They are the same under integration. However diagram 
(e) contrary to diagram (d) is invariant under the exchange P23. Hence, it has to be weighted with a prefactor 1/2. That is 
why in the explicit formula for Ud one has to subtract from J drgXddfrii) several other terms that are needed to remove the 
overcounting coming from these missing prefactors (see ref. [211] ). 



A(r y ) = (A NN (r ll ) 1 A pp (r tJ ),A NP (r lJ ),A PN (r lJ )) . (30) 

where we have also indicated the projectors used in eq. 1131 

Our nodal equations can be written in a compact way by exploiting the convolution formalism. We make the 
definition: 

C a (ry) - (A(r ifc ) I B(r fci )) a , (31) 

where the subscript a denotes the exchange nature of the node k and, consequently, that of the related vertex 
correction. It can be of type d, e or c. The four components of C are given by 



0«) 


= E f*& 


[ dr k A Nb (r ik )B m (r kj ), 




b=N,P 


In 




= E 


[ dr k A pb (r lk )B bN (r kj ), 




b=N,P 


In 




= E 


f dr k A Nb (r ik )B hp (r kj ), 




b=N,P 


In 


<£ p (*«) 


= E p>>£ 


f dv k A Pb {v lk )B bP {v k] ). 




b=N,P 


In 



(32) 



1. Nodal diagrams 

Using the above convolution formalism the chain equations for N^d, Nd e , N e( j and N ee can be recast in a more 
compact way, as follows: 

N dd = {X dd I N dd + X dd ) d + {X de I N dd + X dd ) e + (X dd I N ed + X ed ) e 

N de = {X dd I N de + X de ) d + {X de I N de + X de ) e + {X dd I N ee + X ee ) e 

N ed = {X ed I N dd + X dd ) d + (X ee I N dd + X dd ) e + {X ed I N ed + X ed ) e 

N ee = (X ed I N de + X de ) d + (X ee I N de + X de ) e + (X ed I N ee + X ee ) e (33) 
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The chain equation for N cc can be written in convolution notation as: 

N cc (r 12 )= /x cc |X cc + N cc -ij +(~\X cc + p) + (~ | -i+N cc -P ) . CU) 



where P is given by: 



P(r 12 )= ^X cc |X cc + N cc -ij . (35) 

These last two equations cannot, however, be applied blindly. According to the RF definitions, cc diagrams have 
exchange-line paths connecting 1 directly to 2 and so, since exchange correlations cannot flip spins, only PP and NN 
components are allowed, whereas the components NP and PN vanish for N cc , X CC ,E CC and P. The degeneracy factor 
(d) here is equal to 2, accounting only for the two possible spin states ( the isospin states of a pair are singled out in 
our treatment ). 

2. Composite diagrams 

The equations for the composite functions are a straightforward generalization of the corresponding equations in 
standard PB-FHNC theory, and are given by 

*Sk(ria) - <&(r 12 ) - A^(r 12 ) - 1 , 

X% e (r 12 ) = X: d (r 21 ) = <&(r 12 )[JV&(r 13 ) + E a de (v 12 )] - JV£(r ia ) , 
Keini) = 92 d (ri2){N: e (v 12 ) + E: e (v 12 ) + [N^(v 12 ) + E2 e (r 12 )} 2 

- d[N^i2) ~ ~4*(r X2 ) + E« c (r 12 )) 2 } - ATf e (r 12 ) , 
Kc(ri2) = ff2,(r ia )[A^(r ia ) - \l{v 12 ) + E«(r 12 )] - N? c (r 12 ) + ±* a (r 19 ) , (36) 
where the g d d distribution function is defined as 

S&M = r 2 (r 12 )exp[iV^(r 12 ) + ^(r 12 )] . (37) 
It is useful to define also the following four distribution functions 

92e(ri2) = gf d (r2i) = N2 e (ri2) + X2 e (r 12 ) , 
(rf e (ri 2 ) = ^ d (r 12 ){iV« (r 12 ) + , 

{gexchfeM = -d 9d < d (r 12 )[N«(r 12 ) -^4(ri 2 ) + E« c (r 12 )} 2 + g2 d (E exch )« e (r 12 ) , 

ffc Q c (r 12 ) = iV c Q c (r 12 ) + X«(r 12 ) - ~£ a (r 12 ) , , (38) 
where the components a = NP and a = PN of (g eX ch)ee are identically zero. 

3. Vertex corrections 

In deriving the vertex corrections we need to distinguish whether the vertex to be corrected is a neutron or a proton. 
Extending the derivation of ref.[5l| we get the following expressions. 



U? = 



E P»[ dr 12 {^[X^(r 12 )-^(r 12 )-^(r 12 )( 5 ^(r 12 )-l)] 

b=N,P ^ n 

&[Xg>(r 12 ) - £f (r 12 ) - SZ b (r 12 )gg\r 12 ) - S^(r ia )(s£?(r ia ) -!)]}+ E n d , (39) 
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for the vertex correction of type d, and 

b=N,P Jn 

- £«>i 2 )gf> 12 ) + S^{v l2 )g^(v l2 )\ 

- C[SZ b (r 12 )(9^ d b (r 12 ) 1) + SZ b (r 12 ) 9 Z b (r 12 ) + ^ e b (r 12 )gf/(r 12 ) + S^> 12 ) ff f> 12 )]} 

- dp N £ J dr 12 {s» N (r 12 )g£ N (r 12 ) + e N (v 12 )W% N (r 12 ) - ~4v(r ia )]} + £f , (40) 

for the vertex correction of type e, where A/" is a shorthand for 

K N c N (r 12 ) = O12) - ^ PN dr 32 X™(r 13 )[X™(r 32 ) + N™ (r 32 ) - \i N {r 32 )\ , (41) 
and the quantity S is defined as 

S ij = ^(N ij + E ij ) . (42) 
Obviously , is given by simply interchanging each superscript p <-> n in the above. 



4- Tuio 6ody distribution function 



The above equations can be solved iteratively. At convergence, the solutions can be used to compute the scalar 
two-body distribution function, given by: 



5c(r 12 ) = £ <?f (r 12 ) 

" a,b=N,P 

+ Cee{(9 dl Xl{v l2 ) + {g eX ch)t e (v l2 ))} . (43) 

In the case of symmetric nuclear matter (N = Z) and of a state independent correlation operator (/|| = /j_), the 
quantity g c (i"i 2 ) recovers the two-body distribution function g(i"i 2 ) of Jastrow theory given in eq. (|24|) . This can be 
easily understood by taking into account (i) the sum rule £<j = 1 and (ii) that £ e g d e and £^g e e of our vertex corrected 
PB-FHNC theory coincide with the corresponding distribution functions g de and g ee of the standard one. 



A. Potential energy expectation value 



The expectation value of a two-body potential of the 1*4 type on the trial function Fjl |0] is given by: 

= ~ J dr 12 {v c (r 12 )g c {r 12 ) + 3i> CT (r 12 )p CT (ri 2 ) + v T (r 12 )g T (r 12 ) + 3v aT (r 12 )g aT (r 12 )} . (44) 

The first term on the r.h.s. corresponds to the expectation value of the scalar component of the V4 potential. We 
discuss, in the following, the remaining three terms. Recall that we are dealing with a polarized system with respect 
to isospin (N ^ Z), but with a strictly non-polarized one with respect to spin. 



1. o\ ■ 02 term 



The correlation operator Fjl has no spin dependence. Since < a\ ■ a 2 > = in spin symmetrical matter, in the 
calculation of the expectation value of u (T (r 12 )tTi • er 2 , the direct terms of the distribution function do not contribute. 
On the contrary, the exchange terms carry the spin-exchange operator and one has to take care of the following spin 
algebra 
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FIG. 6: Left panel: Diagram showing how the new S correlation results from the isospin-state flip of a particle. Right panel: 
An exchange loop passing through two different external points. 



In conclusion we have: 



< (<Ji ■ <J 2 )P a (l, 2) > = <{o x . <r 2 ) 1 + ^ ' g2 > 

= <*^P>=f, (45) 

r(r u ) = ^ {{£pN?{g eX cK)l N {vx2) + (£ P p f(g exch ) p J(v l2 )) . (46) 



2. T\ ■ T2 term 



The operator n -t 2 carried by the r-component of the 7)4 potential requires a specific and new PB-FHNC treatment, 
when dealing with a correlation operator of the Fjl type and a N 7^ Z matter. We begin by calculating the isospin 
matrix elements. The direct terms are: 

(NN I v T (n ■ t 2 ) I nNN) = (PP \ v r {n ■ r 2 ) | PP) - v T (r 12 ) , 
(NP I v T (n ■ r 2 ) I NP) - (PN I v T ( n ■ r 2 ) | PN) = -v T (r 12 ) . (47) 

and the exchange terms are: 

(NN I v t (ti ■ t 2 )P t I nNN) = (PP \ v T (n ■ r 2 )P t au \ PP) = v T (r 12 ) , 
(NP I v T (n ■ t 2 )P t I NP) = (PN I v T (n ■ t 2 )P t I PN) = 2v T (r 12 ) . (48) 

The second row of the above equation deserves particular attention. Due to the fact that we have e.g. | NP) in the 
ket and (PN \ in the bra, we have: 

1. a different kind of correlation reaching the external point 1 or the external point 2. Under the assumption of 
eq. (H31) there is onyl one type for such new correlations, and we denote it as (5-correlation: 



hsinj) = /||(ry)/x(rii)-l. (49) 

2. the exchange loop passing through 1 and 2 is made of two cyclic nodal functions one of type P and the other of 
type N. 

It follows that the g-distribution for these matrix elements, which we denote by gs{ri 2 ), has to be built by solving 
appropriate PB-FHNC equations, taking into account the above two properties. We have: 

9rM = lL(r 12 )+ £ {pl[(0 2 9rAri2)+QC(g a d a e(ri2)+g a eS(ri2)) 

P \ a=N,P 

+ (C) 2 ((g^r)Te(ri2) + (ge X c h re(ri2))}} (50) 



- E p*p* [ti$gZ(ri2) + a e a 5f(ri 2 ) + C«(ri2) + CC(gMr)Z(v 12 )) 



=N,P 
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where a labels the isospin conjugate of a, namely \ a) — t x \ a) and the mixed distribution function gs is given by 

1 



gs(r 12 ) = -4dp N p P £ eS fl 



2.f2 JV™/>12) 



(51) 



In the case of N = Z matter and /ii = /x, one has p^ = pp, £, N = X P with the result that gg = 2g a and 
consequently that g T = 3g a . 

The nodal function Ngg P (r^) is defined by the convolution: 

ATJF>12) = (X 4d | N dS + X dS ) d + {X Se | N dS + X d s) e + (Xs d | N e ,5 + X e s) e (52) 

where the nodal and composite vector functions of the type dS and eS have only the two components which specify 
the isospin state of the particle related to the external label d or e. The two chain equations are: 

N 5d (r 12 ) = ((X M | (N dd + (X dd ) d + ((X 5e | N dd + X dd ) e + (X Sd \ N ed + X ed ) e 
(N 5e (r 12 ) = (X Sd | (N de + (X de ) d + ((Xie | (N de + (X de ) e + ((X 5d | (N ee + (X ee ) e , (53) 

and the composite functions are: 

X b Sd (r 12 ) = 9 b sd (r 12 ) - N% d (r 12 ) - 1 , 

X| e (r 12 ) = X b 5 (r 21 ) = g b d (r 12 )[N b e (v 12 ) + E b Se (r 12 )} - N b Se (v 12 ) , (54) 
with b = N,P and gs d given by 

g b Sd (v 12 ) = /||(r 12 )/ ± (r 12 )e( JV «( r -)+ x ~^)). (55) 

The two components Nf Scc , with a = NN, PP , entering eq. (f5"Tj) . can be given in terms of the following convolution 
equation 

N«5«5 cc (ri 2 ) = (^Xscc | X d - cc + Nscc - -j^J ^ + ^ | Xscc + Ps 

~|~+N«cc-P«) , (56) 
where the two-component vector functions N,5,5 CC , Xa,s cc and Pa are given by 

N 5cc (ri 2 ) = (x Scc | X cc + N cc - + f-~ | X c 

*&c(*ia) - <7&(ri 2 )[iVJ* cc (r 12 ) - ^ ( ria ) + £& c (r 12 )] - JV£ c (r 12 ) + ^ a (r 12 ) , 
P«(r 12 ) = ^X cc | X cc<5 + N CC(5 - . (57) 

In order to complete our PB-FHNC set, we need to define Ug so as to define £g e : 

6e = , (58) 



Us = 



£ P" I dr 12 {e d [X b 5d (r 12 ) - E b d (r 12 ) - S b Sd (r 12 )(g b d (r 12 ) 1)] 

u—N.P 

+ e e [X b s e (r 12 ) - E b g e (r 12 ) - S b g d {v l2 )g b e {v 12 ) - S b Se (r 12 )(g b d (r 12 ) - 1)]} + Eg , (59) 



where: 



9 b Se (r 12 ) = N b Se (r 12 ) + X b e (r 12 ) . 

S b t (r 12 ) = \(N b t (r 12 ) + E b t (r 12 )) (t = d,e), (60) 
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3. (<7i ■ (T2) (n ■ T2) ienn 

The distribution function can be easily calculated by using the expressions derived in the previoues two sub- 
sections. 

9,r(ru) = ± {pU€) 2 (9e*c h )Z N (ri2) + Pp(e) 2 (9e^) e T(r 12 ) + gs) ■ (61) 

B. Kinetic energy expectation value 

In this section we calculate the exp ectation value of the kinetic energy (K) /A using the Jackson- Feenberg identity, 
following the procedure shown in [261 ] and further discussed in [2l|. It is given as a sum of three terms: a term giving 
the fermi energy (Ep), a term accounting for two-body contributions (E2) and a term for three-body contributions 
(E3). The fermi energy is easily expressed as: 



Ep — — 



2 U 2 1 fc2p 



1 ^ h 2 k 2 N 1 ^ h 2 k 



N ^ 2m Z ^ 2m 

filled N states filled Z states 

The two body contribution to the kinetic energy can be split into two parts: the first one accounting for the contribution 
coming from V 2 acting on the correlations f\\ and f± and the second one accounting for exchanges. The resulting 
expression for the case of the Fjl ansatz and a N ^ Z matter is: 



E2 — E% + Ef , 

E 2 F = -^-p[ dr[ 5 f>)+ 3 f>)]v 2 l n/|| (r) 

- ^- p Jjr[ g r(r)+ 9 r(r)] V 2 ln/ ± (r), 
Ef = jf«fr(fl5^(r)-l) (jv% 2 (r)-« W Wv 2 ^(r) + 2e^V%(r)) 

~ ^7 L ^ " 1} (^^ p2(r) " ^C F «V 2 ^(r) + aef |v 2 £ F (r)) . (63) 

The three body contribution comes from (v£i2 • VI43) and is given by: 



E;l 



\- Pa 



imd ^ 

a=N,P 



dr / dr'v£ a {r) ■ v£ a (v') 
P Ja Jn 



x (^(r)-l)(^(rO-l)e e ° 



N a cc {v - r') - -4(r - r') 



9%{* 



(64) 



C. Single particle excitation spectrum 

The single particle potential of nuclear matter is calculated by applying the method devised in ref. [27j to PB- 
FHNC theory and N 7^ Z matter. This will allow us to evaluate the single-particle neutron and proton potentials in 
neutron rich matter. 

It is convenient to calculate, as in ref. [13] , the particle-hole excitation energy rather then directly the single particle 
excitation, mainly because the | ph) state has the same number of particles as the ground state. Let us consider 



e a (p,h) = (( Po h Q I H I Pah a ))p a>ha - E , (65) 

where (ph | H | ph) is defined in eq. ©, the label a specifies the isospin nature of the excitation, and () stands for 
average over the directions of p and h. As discussed in section UlI Al the variational estimates of diagonal states are 
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maintained after orthogonalization. The single particle excitation e a (q) can be obtained from e a (q, k F ) = zte a (q)^fe Fa , 
where the upper sign is for particle state {q > kp) and the lower one for hole states (q < kp), and 

E p a dE 
£Fa = ~A + A^i> 

= eo + Pa^p ■ (66) 

The single particle spectrum is related to the real part of the nuclear optical potential by 

e a (q) = ^ + U *(<1)- (67) 

One can eliminate q from e a (q) and U a (q) to obtain an energy-dependent U a (e). Perturbative corrections to e a (q) 
and U a (q) include coupling with two-particle one-hole states | p',p",h) for q > kp or two-hole one-particle states 
| h',h",p), for q < kp, giving a width to | q). 

The particle-hole state | p a h a ) is generated in PB-FHNC theory by introducing the following density matrices 

la(rij;qa;k Fa ) = l a (rij) + /S. a (vij;q a ; k Fa ) , 
A a (r y ; q a ;k Fa ) = ±— ]T ^^(j) T —- £ K(^)Mj) , (68) 

q ° neq a kFa nGk J?Q 

where the summations are extended to the iV q states of the shell T?q\l (2rn) and the iVk F states of the shell 
tfklJ{2m). 

The cluster diagrams contributing to e a (q, h) can be obtained by substituting in all the allowed diagrams of Eq one 
and only one Z a -line with a A a -line, for all the Z a -lines of the diagram. To sum up the resulting series of cluster terms 
one can use the following algorithm 

1. modify the density matrices in the following way 

l a -> l a + xA a , (69) 

where x is a smallness parameter and serves to take only one A-line at time in each diagram; 

2. solve the PB-FHNC equations with the modified density matrices; 

3. compute the energy expectation value e a (p^, pp;q;kp;x), as for the ground state energy, with the Fermi energy 

given by 



h 2 

e Fa (q,kp)^—(±q 2 Tk Fa ); (70) 
Ira 



4. compute the particle -hole excitation from 



e a {q,k F ) = —e a (p N ,p P ;q;kp;x)\ x=0 . (71) 

Note that the discrete character of PB-FHNC implies that e Fa will lie in between two energy shells, that we denote 
as e Fa - and e Fa +- It follows that the particle energies will be extracted from e a (p a ,fcFa-) and the hole states from 
e a (h a , k Fa+ ). Therefore, one also need to compute e a (k Fa+ , k Fa _). 

The neutron and proton effective masses are given by the derivatives 

= j = 9U a (e) 
m de 

Enhancements of m*(e) will correspond to flattening of U a (t) around e ~ e Fa which, most likely, will happen only 
after having added the perturbative corrections [28j. 
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Potential Approx. 


A 


Ef ree 


PE 


KE 


E 




Jastrow 


2060 


22.136 


-43.595 


28.455 


-15.150 


AT4' 


F.JL 


2060 


22.136 


-44.053 


28.016 


-16.090 




F 2 


oo 


22.107 


-44.163 


28.208 


-15.955 




F 4 


00 


22.107 


-44.756 


28.587 


-16.169 




Jastrow 


2060 


22.136 


-27.545 


29.499 


1.954 


(AV8') 4 


F.JL 


2060 


22.136 


-30.070 


31.722 


1.652 




F 2 


oc 


22.108 


-28.569 


30.152 


1.583 




F 4 


oo 


22.108 


-31.146 


32.055 


0.909 



TABLE I: Results for symmetrical nuclear matter at p = 0.16 for different correlated models. The first two rows for each 
potential are obtained for 2060 nucleons in a periodic box by using the vertex corrected PB-FHNC equations. We have used a 
grid with 60 points in each direction. The third and fourth rows report the results obtained with FHNC/SOC equations in the 
thermodynamic limit for the F2 and F 4 models. The energies are in MeV 



A 




PE 


KE 


E 




28 


22.427 


-43.994 


28.311 


-15.682 


-0.888 


76 


21.231 


-44.890 


27.119 


-17.771 


-0.938 


108 


21.277 


-44.933 


27.163 


-17.769 


-0.932 


132 


21.996 


-44.304 


27.895 


-16.409 


-0.935 


2060 


22.136 


-44.053 


28.016 


-16.090 


-0.940 



TABLE II: PB-FHNC Results for the AT4' potential at po- For the smaller systems, 26 neighbor cells were summed over; our 
grid had 60 points in each direction. A_b gives the difference in energy with respect to the Jastrow case. The energies are in 
MeV 



VII. RESULTS 

In this section we present and discuss vertex corrected PB-FHNC calculations, performed with the Fjl model, 
under the parallel-antiparallel approximation of eq. (fT3|) for the ATA' and (AV&'\ potentials. 

We encoded our PB-FHNC scheme as an extension of the code already used in [18| and we refer the reader to that 
paper for details of the numerical techniques used. Whenever possible, we made use of standard libraries (e.g. fftw3) 
and routines documented elsewhere (e.g. the ODE integration routines from [111). We tested our double-precision 
code with different compilers and different optimizations, always obtaining consistent results. 

A. Comparison of various correlated models 

We tested our vertex corrected PB-FHNC scheme by comparing the results for SNM with the simple Jastrow ansatz 
with those of ref [l8[ . The results obtained for the expectation values of the kinetic energy, KE, the potential energy, 
PE and the total energy, E, and displayed in the first and the fifth row of Table U coincide within five digits with those 
of ref [lH . Note that this check is not at all trivial because it follows from the fulfillment of the sum rule ^ = 1 and 
of the relations £ f ..qu P , = g^ e and Qg ee = 9ee > where g^ e and g^ e are the not vertex corrected PB-FHNC distribution 
functions of ref. fl8j |. 

We have analyzed the quality of our proposed variational model Fjl by comparing it against the results obtained 
with the F2 and F4 models for SNM at po ■ As shown in Table U the longitudinal isospin-dependent model correlation 
Fjl improves considerably the Jastrow ansatz for the AT 4! potential. In addition, it gives equally good energy results 
as those of the F2 model. The effectiveness of r 2 (l)r z (2) correlations as compared with the t\ ■ T2 ones is confirmed 
by the results obtained with the (AI/8')4 potential which has a stronger a ■ a dependence. In fig we show our 
correlation functions computed at p = 0.16. As expected from the Pauli exclusion principle, at r = there is a 
stronger antiparallel correlation and weaker parallel one. 

Tables [TT1 and Mil show the dependence of the energy results on the number of nucleons in the periodic box. One can 
see that with A = 132 finite size effects are still large. They become totally negligible at A — 2060 (see also ref.|18||) 
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FIG. 7: AT4' correlation functions at p = 0.16. The parallel and antiparallel correlations of the Fjl model are compared with 
Jastrow correlation. 



A 


Ef ree 


PE KE 


E 


A B 


28 


22.427 


-30.937 31.971 


1.034 


-0.270 


76 


21.231 


-30.618 30.782 


0.164 


-0.250 


108 


21.277 


-30.590 30.826 


0.235 


-0.249 


132 


21.996 


-30.187 31.559 


1.372 


-0.256 


2060 


22.136 


-30.070 31.722 


1.652 


-0.302 



TABLE III: Results for the (AV%')i potential at po- See captions of table ITT1 



B. Clustering at subnuclear densities 



Although our variational ansatz does not allow an explicit clustering of the nucleonic matter, there are strong 
indications for such clustering phenomena for density below O.lpo- To analyze these indications we have studied the 
behavior of total energy and of the pair distribution functions of SNM at p — 0.1p hi two different regions of the 
healing distance d. At small d (d ~ 1 fm) the correlation functions /|| and f± are below 1 for r < d, and the system 
does not show any clustering phenomena. On the contrary, as shown in Fig [5J at large d (3 fm < d < 3.35 fm) 
the correlation functions have pronounced peaks, occurring roughly at the same value of the interparticle distance, 
irrespective of the value of d, as typically happens in clustering phenomena. One can see from Table IIVI that the 
variational energy gets lower in the region of large d reaching a minimum around d ~ 3.3 fm. 

In table IIVI we show the energy of a system of 2060 nucleons (SNM) computed with different healing distances. For 
all these calculations we enforced the same level of convergence of our FHNC equations. For d > 3.35 fm we were not 
able to make our equations converge any longer. For comparison we show the energy obtained with d — 1 f m which 
is the largest healing length that does not produce the clustering effect. 

In figure [9] we show our results, obtained with d — 3.30fm for the various g components and for gdd, 5deandg 00 . It 
is easy to see that indeed in the channel np there is a clear evidence of clustering occuring. 
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FIG. 8: AT4' correlation functions at p — 0.016 for different healing distances (d) either for /y (marked with p) or fx (marked 
with a). We show for comparison two correlation functions whose healing distance was set to 1 so to keep their value everywhere 
below 1. 



dfm 


1 


3 


3.1 


3.2 


3.3 


3.35 


EMeV 


-0.73 


-2.634 


-2.787 


-2.950 


-3.087 


-3.097 



TABLE IV: SNM energy per particle computed at p = 0.016 using the AT4' potential. Our grid was set up to have 45 points 
in every direction. 



Our results strongly indicate that at such low densities variational functions allowing for a nucleus embedded in a 
neutron fluid would have a lower energy with respect to those describing an homogeneous fluid. In the following we 
have always forced the system to behave as an homogeneous fluid. 



C. Symmetry energy at subnuclear densities 



One of the main advantages of a t z \t z % form of the isospin correlations is that one can easily compute the energy 
expectation value of nuclear matter with N Z, provided that both N and Z are magic numbers (i.e. they correspond 
to shell closure). For this reason, we cannot keep A = N + Z fixed. The average value of A which we find more 
convenient because it is sufficiently large to reduce finite size effects and allows for a quite large number of admixtures 
with the smallest fluctuations (AA) is the magic number 1898. In table [V] we report the number of protons and 
neutrons for each admixture, the percentage Z / A = 0.5(1 — a) and the energy results of our calculations - performed 
using AT4' potential - at 4 different densities. 

We used our results to check whether the standard way of fitting asymmetric admixture energies (i.e. using a 
simple quadratic fit) can be reliably used at sub-nuclear densities. We used Mathematica to fit our results using a 
6th degree polynomial as a prior. As expected we get null coefficients for the odd power terms. We also get non 
zero coefficients for the 4th power term. Such coefficients however are always negligible, being, at most, one order of 
magnitude smaller than the 2nd power term one. 
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FIG. 9: AT4' gfunctions at p = 0.016 for d=3.30. 
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Q 2 


E(W- 2 p ) 
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E(p = 0.: 





1898 


1898 





0,00% 


1,000 


1.26 


4.94 


25.55 


57.13 


54 


1850 


1904 


6 


2,84% 


0,890 


1.16 


4.28 


20.78 


47.52 


114 


1790 


1904 


6 


5,99% 


0,775 


1.05 


3.62 


15.95 


39.30 


186 


1694 


1880 


-18 


9,89% 


0,643 


0.93 


2.86 


10.44 


29.00 


294 


1598 


1892 


-6 


15,54% 


0,475 


0.78 


1.91 


3.44 


15.84 


406 


1502 


1908 


10 


21,28% 


0,330 


0.65 


1.08 


-2.64 


4.40 


514 


1382 


1896 


-2 


27,11% 


0,210 


0.55 


0.42 


-7.55 


-4.88 


730 


1174 


1904 


6 


38,34% 


0,054 


0.42 


-0.43 


-13.84 


-16.79 


874 


1030 


1904 


6 


45,90% 


0,007 


0.38 


-0.70 


-15.88 


-20.67 



TABLE V: Number of protons and neutrons for different admixtures that we studied and results at different densities. 



In figure [TU] we show our results and our fits for the energies (E) of admixtures with different a 2 . In table I VII we 
give our best estimate for the symmetry energy (S) at different densities computed using the AT4' potential and a 
purely quadratic fit. The value obtained at po is somewhat larger than the experimental results of S ~ 36MeV; this 
is expected due to the phenomenological nature of our potential. 
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FIG. 10: Energies per particle (E/A) of admixtures with different a 2 at different densities. 



p(fm" 3 ) 


10" 2 po 10"Vo Po 0.25 


S (MeV) 


0.89 5.66 41.62 77.96 



TABLE VI: Symmetry energy computed at different densities using the AT4' potential. 



D. Equations of state for AT4' 

In figure [IT] we show our results for the energy of SNM computed at different densities using the AT4' potential. 
We find that we can nicely fit our data by: 



Esnm(p) = E q + a(p - Peq ) 2 + b(p - p eq fe^-P^ (73) 

where E = -22.51MeV, p eq = 0.33 far 3 , a = 220MeVfm 6 , b = -1.56MeVfa 9 and 7 = -5.570fm 3 . The EOS 
for asymmetric matter can then be written, fitting the data for different a, as: 



E(p, a) = Esnm(p) + C s (-?-) a 2 (74) 

\PeqJ 

where C s = 82.86 MeV and 7s = 0.913. 

In figure [T2l we show our EOS at subnuclcar density computed in two cases: PNM and 10% protons. We show in the 
same plot the BPS EOS [30| as a useful comparison. Quite comfortingly our results show quite a good agreement at 
the edge of the inner crust. This agreement is obviously not preserved at lower densities due to the absence of clusters. 
In the same diagram we also show a point computed subtracting, from the energy of the pure gas, the binding energy 
of the corresponding nucleus (as reported by |3l|) obtained using the semiempirical mass formula. We used different 
polynomials to fit E(p) and hence to derive P. The errorbars show our best estimate obtained using a cubic spline 
interpolation. 



E. Particle-hole 



In fig[T3] we show our results for the single particle energy in different test cases as a function of q = p — kf. The 
values at q = are given by the corresponding chemical potentials e^a given in eq. [SH The correlation effects can 
be viewed by comparing the F4 results with the corresponding Fermi gas estimates. Our results at po and a = are 
in reasonably good agreement with the FHNC/SOC calculations of ref [28| obtained with the Urbana V14 + TNR 
interaction. One can see that at lower densities the effects on the optical potential due to the asymmetry are much 
reduced. 
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FIG. 11: Energy per particle of SNM computed at different densities using the AT4' potential. For low densities we forced the 
system to be homogeneous. 
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FIG. 12: Energies (E) of admixtures with different a 2 at different densities using the AT4' potential. 

VIII. CONCLUSIONS AND PERSPECTIVES 

We have made a first step towards the development of a technique, based on CBF theory, allowing for a study of 
the NS crust from first principles. In particular we have developed a theoretical framework suitable for studying NS 
crustal cells using V4 potentials and we have applied it, using the AT4' potential as a test case. 

Our results are promising in several respects. First of all it is shown that using only the third component of the 
isospin dependent correlation is a good enough approximation. We are in a good position, from the variational point 
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FIG. 13: Single particle energy e(p) for N — Z and N = 2Z at p = po (upper panel) and p — O.lpo (lower panel) computed 
using the AT4' potential. The figures report also the fermi gas estimates. Energies are in MeV and momenta in fm _1 
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of view, to insert a nucleus in our system as discussed in the introduction. We are also in a good position to rewrite 
the FHNC/SOC scheme, and the CBF perturbative corrections, using this simplified operator to treat the isospin. 
This would enable us to use more realistic potentials (particularly those with tensor interaction) and hence to refine 
our results. Notice that standard FHNC/SOC does not allow for asymmetric matter. 

Moreover the vertex corrected theory developed here can be readily used to include superfluid effects, thus improving 
the accuracy of our description of the crustal matter. 

Finally the extension of PBFHNC theory to the treatment of asymmetric matter is essential to deal with the crustal 
matter, thus, enabling a fully self-consistent description of the neutron star equation of state. 

This work was partailly supported by CompStar, a Research Networking Programme of the European Science 
Foundation. This work was partially funded under MIUR PRIN fund "Fermi Hypernetted Chain and Quantum 
Monte Carlo studies of nuclei and nuclear matter with applications to the astrophysics of neutron stars" , and National 
Science Foundation grand PHY0757703. 



Appendix A: Computation of correlations 



To compute the correlation function needed for our calculation, we have to solve a set of differential equations which 
can be derived by minimizing the expectation value of the energy given by the lowest order diagrams, as shown in 
[25| . The expectation value of the energy is: 



1 



t 



(v c + vr) [p^ 1 - -iff + Pp 1 - ~e P - 3 k + v aT ) - Pp i P + -put 



-N 



PpPn 
P 



I dr{^( V .^) 2 + /] 
Jn rn 



}• 



(Al) 



Minimizing this expression with respect to /m and f± gives: 

-£(v/l) 3 + /l[i>c-'. 



T - £ P £ n N {V T + iv aT )] = , 

(V/||) 2 G + v/|| (VG)1 + f\\ [(v c + vr) G - 3 (iv + v aT ) (p% + p% — G)] = . 



where: 



G = pjf 1 



'-£ 



N 



PP 1 



(A2) 



(A3) 



Without any loss of generality, we can divide each term of the second equation by p 2 and redefine G accordingly. In 
the following, we set C, a — p a / P ■ 

These two differential equations need to be solved bearing in mind that / should heal smoothly to one at some 
distance d. For easily solving them, it is convenient to start from r — after having redefined the variables: 

4>± = rf± , 
^11 = rVG/, 



with: 



0x(d) = d 4>' ± (d) = 1 . 



1 dG'(d) 



i&„(d) = d^/G{d) ^(d) = VWj+z v ^ 7T 

These two new functions are defined so that their value at the origin is zero. Since there are then three boundary 
conditions, we need to introduce two lagrange multipliers A to ensure that they are all satisfied simultaneously. This 
leads to: 



-^y +^ll 



^0" + 0J_ [V c - V T - £ p £ n (V T + 3U CTT )] = \4,(f>± , 

h 2 ( 1 G" I G' 1 G l2 \ 
m \ 2 ~G~ ' rG 4 ~G T J 



(v a +v T )G-3(v a +v^)(C+< 2 p ~G) ^ ^ 
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(A4) 
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We solve these two equations using a standard adaptive-stepsize Bulirsch-Stoer method, varying d so as to minimize 
the energy and iterating to evaluate the As. The equations have the form 



0" + (a(r) - A) (f> = , 
and we can adjust lambda by adding, at each iteration, a SX defined as: 



SX 



(A5) 



(A6) 



where the quantities denoted with a subscript C are those computed numerically at the previous step, and those with 
a subscript T are theoretically derived boundary conditions which need to be satisfied at r = d. 



Appendix B: Some standard quantum mechanics results 

As a useful reference we recall that: 

A[ipl ■ ■ ■ ipx]0(xi, • • • , x n )A[ipi ■ ■ ■ ipn] 

can be rewritten as: 

1 N 

-j= X] ' ' "^n]0(xx, ■ ■ ■ ,x n )A[ipi ■ ■ -iPn] ■ 
The anti-symmetrizing operator can be written in terms of the two-particle exchange operator P: 



A = l-J2 P v+ ( P H P ik + PikPkj) + ( P « Pj W + P ^ P oi + p u p jk) - {ijklloop} 

i<3 i<j<k i<j<k<l 



(Bl) 



(B2) 



(B3) 



Applying eq. IB31 it is easy to check that the fermi-gas n-body correlation function is given by 



^ G (n,...,r„) = 1- ^a(r M ) 2 + ^(riMrjMrk,- 

a,(i<j) a,(i<j<k) 

and, in a more compact form, by the following determinant 



(B4) 



9n (ri,...,r„) 



1 3E 4(1,2) JE.4(1,3) ... ^£„4(1, 

3E„4(2,1) 1 3 £.4(2,3) ... i£„4(2 : 

3 £«*«(«,!) i£ a 4(A2) i£ a 4(n,3) ... 1 



where £ a (hj) and £ a (rij) are defined in eq. (|16p . the summations run over the isospin states, (a = N,P) and 
e? = 2. The products of the £ a (hj) operators also implies the matrix elements of the relative spin-isospin states (all 
the spin-isospin states of the particle in a loop must be the same). In infinite matter £ a ( r i,j) reduces to 



40*i,j) = £{x = k Fa rij) = -3 [sin(x) - xcos(a;)] 



(B5) 
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